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Abstract 



& An algorithm for integration of polynomial functions with variable weight is 

i 1 considered. It provides extension of the Gaussian integration, with appro- 

Qh priate scaling of the abscissas and weights. Method is a good alternative to 

usually adopted interval splitting. 
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^ 1. Introduction : the integral 
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In many areas of the physics and chemistry, the following integral emerges: 
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where gr(a;) is a smooth function and < a < oo a real parameter. Here we 
assume that function g(x) can be successfully approximated by the polyno- 
mial. For a fixed value of a the standard method for numerical evaluation of 
such an integral is Gaussian integration or similar closely related algorithm. 
K> The difficulties caused by the integral of the form (fTl) can be show by 

comparison with other similar examples. If we change lower integration limit 
from 1 to zero, one can easily remove parameter a from the algorithm by the 
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substitution x = at: 
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dt. 



Parameter a now appear in the function g, and modified (because of 1 in the 
denominator) Gauss-Laguerre algorithm can be used. 

Similarly, if we remove 1 from the denominator, the substitution z = 
(x — l)/a transform integral into form: 
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easily integrable numerically with standard Gauss-Laguerre algorithm. 

However, when both the lower integration limit is 1, and 1 is present in 
the denominator, parameter a cannot be eliminated from Eq. ([l]): it always 
appears either outside function g in non-linear way, or enters the limit (s) of 
integration. This does not prevent us from calculating abscissas and weights 
of the Gauss-like quadrature for any fixed value of a, say a = 1, a = 3 
or a = 1/4. The general idea of the algorithm presented in this article, is 
to use Gauss-like quadrature with abscissas and weights being functions of 
parameter a. 

Standard method to handle integrals of the form (fTl) is to split integration 
interval into at least two: l<x<£, £ < x < oo. If the value of £ is chosen 
properly as a function of a, then in the first interval we have: 
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that is, weight is nearly linear (or polynomial if using higher order expansion). 
In this interval we can use Gauss-Legendre integration. For x > £, 1 in 
the denominator can be omitted, and Gauss-Laguerre quadrature apply. In 
practice, more interval subdivisions are required [H [21 E] . 

2. Case with fixed parameter 

We start analysis with simplest case of a = 1 in Eq. (pi). The integral 
becomes: 
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For g(x) = x n , where n is non- negative integer, we have founcr] that 
moments are equal to: 
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where Li^ is the polylogarithm. Our goal is to find weights, Wi, and abscissas, 
Xi, of the Gaussian quadrature formula: 
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where P n is polynomial of the order up to 2N gauss — 1. 

Convenient method [4] to find weights and Gauss points uses orthogonal 
polynomials Q n related to given weight: 
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Unfortunately, as noted by [5], for infinite integration interval calculations 
still require arbitrary precision calculations. Problem is increasingly ill- 
conditioned numerically. We have found, that the best method is to solve 
system of equations for unknown polynomial coefficients in eq. Q progres- 
sively. We use already found coefficients for polynomials of the order n — 1, 
and Laguerre polynomial L n coefficients as a guess starting points. System of 
the equations is then solved numerically using arbitrary precision arithmetic. 
Results were verified using method of [Bj. 

Once orthogonal polynomials are found, abscissas are zeros of the Q n , 
and weights are equal to [5]: 
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1 Although I do not have a proof of this formula, it was verified using Mathematica up 
to n = 23. 



Another equivalent formula for weights is [3]: 
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Knowledge of abscissas and weights for integral with a = 1 is crucial, 
because of the approximate scaling found and used in the next section to 
derive more general results for o^l. 

3. Scaling of abscissas and weights 

For integral pj) with function g{x) being equal to x n we can provide 
formula similar to (J3J), generalized to case a^l: 
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Knowledge of the moments facilitates calculations of abscissas and weights, 
but in practice polylogarithms present in eq. (J5l) often is not directly avail- 
able. One must calculate moments by means of direct numerical integration. 
In principle, using (J5J), one can find orthogonal polynomials (and Gaussian 
quadrature abscissas and weights as well) analytically. However, formulae be- 
come ridiculously complicated already for n > 2, and we restrict discussion 
to case of n — 1. 

Unfortunately, we were unable to derive general formulae for orthogonal 
polynomials Q n or find coefficients of the three-term recurrence formula. 

3.1. Single Gaussian point quadrature 

While the case of single Gaussian point is not interesting from practical 
point of view, it gives some insights into planned procedure due to very simple 
formulae. Quadrature is: 
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where: 



W i = aln(l + e- 1/a ) = aLi 1 (-e~ 1/a ), x x = 1 + a ^ & u\ (7) 
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Figure 1: Weight (left) and abscissa (right) and for single point Gaussian quadrature as 
a functions of a, and its asymptotes (dashed). 



Integration algorithm based on abscissa and weight provided above is exact 
only for constant functions. Weight and abscissa are shown in Fig. [TJ Func- 
tions Lifc(— e~ 1//a ) present in the formulas above are somewhat pathological. 
All derivatives vanish at a = 0: 
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Therefore, the function ip{a) defined as follows: 



ip(a) 
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is an example of the infinitely differentiable function which is non-analytic at 
a = 0. Therefore, it cannot be approximated by the polynomials near a = 0. 
For a-loowe have: 
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where ( is Riemann-zeta function. Both functions W\{a) and Xi(a) approach 
their asymptotes very slowly, cf. Fig. [TJ 

3.2. General Gaussian quadrature 

For o^l, due to extremely complicated formulae, discussion will be lim- 
ited to numerical results. It is observed, that functions Wi(a) and X^a) be- 
have like functions Wi(a) and Xi(a) discussed in previous subsection. Rough 



approximation for abscissas and weights (arbitrary order) can be combined 
from results for a — 1, and single point quadrature as follows: 

Xi(o)~l + oAi, l^(a)~(5 l (-l + aln(l + e 1/a )) (10) 

where A* = 1 + X^l) and S t = Wi(l)/(-l + In (1 + e)). Values of Wi(l) and 
-Xj(l) are simply weights and abscissas for a = 1 quadrature, cf. Sect. [2j 

Example results are presented in Fig. [2] using lines. Gauss points behave 
as expected for a — > 0, that is, they concentrate near x — 1. For a — >• 
oo abscissas spread from x = 1 to infinity. Weight also behave correctly, 
scaling linearly with a — > oo and approaching zero as a — > 0. For a — > 



0, formulae (10) are progressively better approximations. Note however, 
that (10) do not provide correct asymptotic behavior for a — > oo. This 
would require knowledge of analytical formulae for Wi(a) and Xi(a), or its 
asymptotic expansion at least. 

Parameterization ( 10 ) provides excellent starting guess pointaj for nu- 
merical calculations of Wi(a) and X,(a). 

Typical numerical result is shown in Fig. [2] using bullets. Apparent ac- 



curacy of the scaling (10) is misleading, because (i) logarithmic scale hide 



errors (ii) very high accuracy is required for Gaussian quadrature to work. 



Therefore quadrature obtained from ( 10 ) is not expected to be very accurate 



except maybe for a <C 1. Nevertheless, overall picture is appealing, and use 



of scaling, possibly not as simple as (10), seems to be move in the promising 



direction. At present, however, we are forced to use interpolation of numeri- 
cally computed values. For values of a outside interpolation domain we can 
use asymptotes for a ^$> 1, and ([8| for a <ti 1. 



4. Performance of the algorithm 

Approximate value of the integral ([!]) is computed from formula: 

J> gauss 

1(a) ~ J2 Wi(a)P n (Xi(a)) (11) 

t=i 



2 Without guess points it is still possible to find Xi(a), Wi(a) starting with known values 
Xi(l),Wi(l), slowly changing a and using previously obtained values in the next step as an 
another guess. However, this procedure is inherently sequential, while using parameterized 
guess one can do all calculations in parallel. Mathematica script of [5] also do not need 
any guess points, but require 500 seconds to find 10-point quadrature for (pi). Our method 
do the same job in 2 to 3 seconds on the same machine. 
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Figure 2: Numerically obtained Gaussian weights and abscissas for integral of the form 
([!]) (points) compared to formulae (10) (lines). 7-point quadrature is shown. 



where, in contrast to Q, W{ and Xi are functions of a, cf. Fig. [2] 

There are two factors determining accuracy of the method: number of 
Gaussian points, N gauss , and accuracy of the functions Wi(a),Xi(a), i = 
1 . . . N gauss . The latter depends on number of sampling points, N samph and 
algorithm used to interpolate between them. Ideally, we would like to have 
many Gauss points, and very accurate (or exact) scaling of them. In practice 
we should know what is better: large N gauss with small N samp i, or vice versa. 
We address these questions in next subsections. 

4-1. Accuracy as a function of Gaussian points N gauss 

In this subsection we keep number of sampling points used to interpolate 
functions Wi(a) and X;(a) fixed. Number of these functions will be var- 
ied. Naively, from definition of the Riemann integral, we expect that more 
sampling points could result in increased accuracy regardless of the algo- 
rithm used, as long as integral is convergent. On the other hand Gaussian 
quadrature of polynomials is exact up to certain order, but only for precisely 
determined abscissas and weights. Algorithm is then expected to fail (in 
the sense of accuracy) for polynomials except at the grid points. For non- 
polynomial functions Gaussian integration is only approximate functional, 
and accuracy of Wi,Xi is possibly less important. 

Now we attempt to test these expectations using two examples of interest: 

gi(x) = x + x 4 — x 5 , polynomial exactly integrable (12a) 

92(x) = x 4 a/1 + x, function arising from relativistic momentum (12b) 

Number of sampling points has been fixed to 50 per decade per function. 
In the range of 0.01 < a < 100 total of 201 points spaced logarithmically 
were used. Third and eight order interpolation were used for abscissas and 
weights, respectively. 

Three cases are compared: N gauss = 3,7 and 11. Relative accuracy ob- 



tained for polynomial (12a) is presented in Fig. pi Differences are barely 
visible, and error is entirely due to variations of the abscissas and weights. 
Normally, for polynomial of the order < 5 we expect to achieve relative error 
of the order of machine epsilon (Fig. [3j thin line near bottom; axis is placed 
at e = 2.2 x 10~ 16 ). Therefore, this figure shows errors caused by the in- 
terpolation of the abscissas and weights. Error estimate reach huge values 
for small values of a. However, in typical application integrals of the form 
(pi) are calculated over wide range of parameters and summed up. Therefore 
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Figure 3: Relative errors versus a for N gau8S = 3 (dotted), 5 (dashed) and 11 (solid) for 
polynomial (12a). 
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Figure 4: Like in Fig. [31 but for non-polynomial function ( 12b I 



overall absolute error is mainly due to regions with large a where the numer- 
ical value of the integral is largest. Interpolation error visible in Fig. [3] might 



be significantly reduced if required, cf. Sect. 4.2 and Figs. M p 



Test case #2 (12b) provides more realistic task for the algorithm. This 
time number of Gaussian points do matter, especially for a > 1. Thin lines 
show relative accuracy achieved exactly at grid pointqj For a > 1, algorithm 
surprisingly provides accuracy nearly identical to normal Gaussian integra- 
tion. It is however clear from Fig. [4J, that accuracy better than ~ 10~ 10 
cannot be achieved, regardless of the value of N gauss . We need more dense 
grid of points used to interpolate functions Wi(a),Xi(a) in (11), or better 
interpolation algorithm, cf. Sect. 4.2 Again, for a< 1 relative accuracy 



is poor, but integral values are very small in this regime due to exponential 
factor in the integrand of ([!]). For small a one can also consider use of (10) 
instead of numerical values, see Fig. [5j red line. 



3 To get this visual effect, we use original grid to draw thin lines, and irrational step to 
draw thick lines. 
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Figure 5: Relative errors versus a caused in the final integral using increased number 
of interpolation points used to restore abscissas and weights for integrand (12a). Red - 
analytical formula, blue - 50/decade, black -200/decade, green - 1000/decade. 



4-2. Accuracy as a function of sampling points N samp i 

In this subsection we keep number of Gaussian points fixed, N gauss = 7. 
Number of sampling points will be increased, so interpolated functions Wi(a),Xi(a) 
will approach their exact values. For polynomials up to the order 2N gauss — 
1 = 13, integration algorithm is expected to achieve machine precision - in 
limiting case at least. For other functions it should be no worse than ordinary 
Gauss-Laguerre-like integration with fixed weight. 

Typical result is presented in Fig. [5] Let us remind, that in the case 
of exact abscissas and weights for polynomial (12a) we expect to achieve 

2.2 x lCT 16 



machine precision. In Fig. |5j e = 2.2 x 10 , where lower axis has been 
placed. Third order interpolation has been used do compute abscissas and 
eight order to compute weights. Analytical formula (10) is shown in red 



for reference. It is not surprise, that using progressively more interpolation 
points we get better results (Fig. [5]). Using 1000 points per decade (Fig. |5j 
green), relative accuracy is at level of 10~ 13 . 
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Figure 6: Relative errors versus a caused in the final integral using various interpolation 
algorithms used to restore abscissas and weights (using 1000 points per decade) for in- 
tegrand (12b). Red - analytical formula, blue - Hermite, black - splines, green - mixed 



orders (see text). 



4-2.1. Influence of the interpolation order and method 

Interpolation goal is to reproduce abscissas and weights from discrete set 
of points as accurately as possible. Besides number of these points, which 
is obvious factor, secondary important issues are: interpolation method, in- 
terpolation order and location of the points. Influence of various factors on 



relative integration error in the case of function ( 12b ) is presented in Fig. 6 



4000 points spaced logarithmically has been used in the range from 0.01 to 
100, i.e. thousand points per decade. Two standard interpolation algorithms: 
Hermite and spline were used with order 1 (linear), 3 and 8. Noteworthy, 
relative error gain obtained from increase of the interpolation order from 3 
to 8 is as high as six orders of magnitude. Further test has shown, that 
gain is almost entirely due to accuracy of weights. Abscissas are successfully 
approximated even using linear interpolation, cf. green lines in Fig. |6j In- 
terpolation order and method is important mainly if a < 1. For a > 1 no 
differences are visible. 
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Possibly, using non-linearly scattered points and renormalized function 
values we would be able to increase accuracy further without increasing num- 
ber of points. This might be important issue in real-world implementation 
of the algorithm. Too large amount of data can result in cache misses of 
modern CPUs, slowing down calculations. 

5. Concluding remarks 

Novel algorithm for evaluation of the improper integrals with a real pa- 
rameter (fTl) in the weight function has been presented. Method is based on 
Gaussian integration, where abscissas and weights are functions of parameter 
a rather than just real numbers. Both analytical and numerical approxima- 
tions of these functions are considered. Algorithm is shown to be robust and 
useful. 

Presented method might be considered to be an example of rendering 
of the multivariate function into number of functions of one variable. The 
latter can be approximated very accurately, and do not consume computer 
memory, in contrast to tabulations of multivariate functions. Therefore the 
fact, that functions W% are not already known in terms of e.g. elementary 
functions is not serious limitation for modern hardware. 

Possible direction of further research are: (i) search for exact scaling for- 
mulae or their approximations (ii) scaling in limiting cases a — > and a — > oo 
(iii) dedicated interpolation formulae (iv) use of scaling with non-Gaussian 
(e.g. trapezoidal) integration rules. Properties of the related orthogonal 
(non-classic) polynomials are also of interest, particularly three-term recur- 
rence formula and behavior of the leading terms. 
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